# Cleaning the Work Environment ####
rm(list=ls(all=T))
gc()

getwd()

library(sjPlot)
library(multcomp)
library(stargazer)

#### Loading the Dataset ####

s3 <- readRDS(file = "SEEBSS - Study 3.rds")

#### Data Preperation ####

names(s3)[names(s3)=="libcon_self"] <- "Self_Conservatism"
table(s3$Self_Conservatism, exclude = NULL)

table(s3$gender, exclude = NULL)
names(s3)[names(s3)=="gender"] <- "sex"
s3$sex <- 1 + s3$sex*-1 # recoding # 0 - female, 1 - male
table(s3$sex, exclude = NULL) 

table(s3$age, exclude = NULL)

table(s3$education, exclude = NULL)

table(s3$lr_self, exclude = NULL)
names(s3)[names(s3)=="lr_self"] <- "self_left_right"
table(s3$self_left_right, exclude = NULL)

table(s3$pid2, exclude = NULL)
table(s3$pid2_7_TEXT, exclude = NULL)

s3$Partisanship <- NA
s3$Partisanship[s3$pid2==1] <- "Justice and Development Party (AKP)"
s3$Partisanship[s3$pid2==2] <- "Republican People's Party (CHP)"
s3$Partisanship[s3$pid2==3] <- "Other Parties on the Left"
s3$Partisanship[s3$pid2==4] <- "Other Parties on the Right"
s3$Partisanship[s3$pid2==5] <- "Other Parties on the Right"
s3$Partisanship[s3$pid2==6] <- "Other Parties on the Right"
s3$Partisanship[s3$pid2_7_TEXT=="Emek partisi" |
                  s3$pid2_7_TEXT=="Emek Partisi" |
                  s3$pid2_7_TEXT== "EMEP (Emeğin Partisi)" |
                  s3$pid2_7_TEXT== "Hdp" |
                  s3$pid2_7_TEXT== "Komunist Parti" |
                  s3$pid2_7_TEXT== "sol" |
                  s3$pid2_7_TEXT== "Sol parti" |
                  s3$pid2_7_TEXT== "Sol Parti" |
                  s3$pid2_7_TEXT== "Sosyalist" |
                  s3$pid2_7_TEXT== "Sosyalist parti" |
                  s3$pid2_7_TEXT== "Sosyalist pati" |
                  s3$pid2_7_TEXT== "TİP" |
                  s3$pid2_7_TEXT== "Tkp" |
                  s3$pid2_7_TEXT== "Türkiye İşçi Partisi" |
                  s3$pid2_7_TEXT== "Türkiye Kominist Partisi" |
                  s3$pid2_7_TEXT== "Vatan" |
                  s3$pid2_7_TEXT== "vatan partisi" |
                  s3$pid2_7_TEXT== "VATAN PARTISI" |
                  s3$pid2_7_TEXT== "VATAN PARTİSİ"] <- "Other Parties on the Left"
s3$Partisanship[s3$pid2_7_TEXT=="Adalet ve Özgürlük" |
                  s3$pid2_7_TEXT=="Bağımsız Türkiye Partisi" |
                  s3$pid2_7_TEXT== "Bbp" |
                  s3$pid2_7_TEXT== "Demokrasi ve Atılım Partisi (DEVA Partisi)" |
                  s3$pid2_7_TEXT== "Deva" |
                  s3$pid2_7_TEXT== "Deva partisi" |
                  s3$pid2_7_TEXT== "Deva partisi'ne karşı da olumluyum ama politikalarının söylemlerinin netleşmesini bekliyorum" |
                  s3$pid2_7_TEXT==  "Gelecek PARTİSİ" |
                  s3$pid2_7_TEXT== "Saadet partisi" |
                  s3$pid2_7_TEXT== "Yeniden Refah partisi" |
                  s3$pid2_7_TEXT=="Yeniden Refah Partisi"] <- "Other Parties on the Right"
s3$Partisanship <- ifelse(is.na(s3$Partisanship),"Bipartisan or Non-Responsive", s3$Partisanship)
table(s3$Partisanship, exclude=NULL)
s3$Partisanship <- factor(s3$Partisanship, ordered = F) 
table(s3$Partisanship, exclude=NULL)

s3$bipartisan_na <- NA
s3$bipartisan_na[s3$Partisanship=="Bipartisan or Non-Responsive"] <- 1
s3$bipartisan_na[s3$Partisanship!="Bipartisan or Non-Responsive"] <- 0
table(s3$bipartisan_na, exclude = NULL)

s3$akp <- NA
s3$akp[s3$Partisanship=="Justice and Development Party (AKP)"] <- 1
s3$akp[s3$Partisanship!="Justice and Development Party (AKP)"] <- 0
table(s3$akp, exclude = NULL)

s3$chp <- NA
s3$chp[s3$Partisanship=="Republican People's Party (CHP)"] <- 1
s3$chp[s3$Partisanship!="Republican People's Party (CHP)"] <- 0
table(s3$chp, exclude = NULL)

s3$other_right <- NA
s3$other_right[s3$Partisanship=="Other Parties on the Right"] <- 1
s3$other_right[s3$Partisanship!="Other Parties on the Right"] <- 0
table(s3$other_right, exclude = NULL)

s3$other_left <- NA
s3$other_left[s3$Partisanship=="Other Parties on the Left"] <- 1
s3$other_left[s3$Partisanship!="Other Parties on the Left"] <- 0
table(s3$other_left, exclude = NULL)

mean(s3$akp, na.rm = T)
s3$beta_akp <- (s3$akp - mean(s3$akp, na.rm = T))
table(s3$beta_akp, exclude = NULL)

mean(s3$chp, na.rm = T)
s3$beta_chp <- (s3$chp - mean(s3$chp, na.rm = T))
table(s3$beta_chp, exclude = NULL)

mean(s3$other_right, na.rm = T)
s3$beta_other_right <- (s3$other_right - mean(s3$other_right, na.rm = T))
table(s3$beta_other_right, exclude = NULL)

mean(s3$other_left, na.rm = T)
s3$beta_other_left <- (s3$other_left - mean(s3$other_left, na.rm = T))
table(s3$beta_other_left, exclude = NULL)

table(s3$GROUP, exclude = NULL)
names(s3)[names(s3)=="GROUP"] <- "Group"

s3$Coronavirus_Prime <- NA
s3$Coronavirus_Prime[s3$Group == "Treatment 1 - Coronavirus Disease"] <- 1
s3$Coronavirus_Prime <- ifelse(is.na(s3$Coronavirus_Prime), 0, s3$Coronavirus_Prime)
table(s3$Coronavirus_Prime, exclude = NULL)

s3$Terror_Prime <- NA
s3$Terror_Prime[s3$Group == "Treatment 3 - Terrorist Attack"] <- 1
s3$Terror_Prime <- ifelse(is.na(s3$Terror_Prime), 0, s3$Terror_Prime)
table(s3$Terror_Prime, exclude = NULL)

s3$Drunk_Driving_Prime <- NA
s3$Drunk_Driving_Prime[s3$Group == "Treatment 2 - Drunk Driving Crash"] <- 1
s3$Drunk_Driving_Prime <- ifelse(is.na(s3$Drunk_Driving_Prime), 0, s3$Drunk_Driving_Prime)
table(s3$Drunk_Driving_Prime, exclude = NULL)

table(s3$treatment_1_1, exclude = NULL)
s3$treatment_1_1 <- 8 - s3$treatment_1_1
s3$treatment_1_1[s3$Group=="Control"] <- 0
s3$coronavirus_prime_fear <- s3$treatment_1_1
table(s3$coronavirus_prime_fear, exclude = NULL)
table(s3$Group, s3$coronavirus_prime_fear, exclude = NULL)

table(s3$treatment_3_1, exclude = NULL)
s3$treatment_3_1 <- 8 - s3$treatment_3_1
s3$treatment_3_1[s3$Group=="Control"] <- 0
s3$terror_prime_fear <- s3$treatment_3_1
table(s3$terror_prime_fear, exclude = NULL)
table(s3$Group, s3$terror_prime_fear, exclude = NULL)

table(s3$treatment_2_1, exclude = NULL)
s3$treatment_2_1 <- 8 - s3$treatment_2_1
s3$treatment_2_1[s3$Group=="Control"] <- 0
s3$drunk_driving_prime_fear <- s3$treatment_2_1
table(s3$drunk_driving_prime_fear, exclude = NULL)
table(s3$Group, s3$drunk_driving_prime_fear, exclude = NULL)

s3$priming_fear <- rowSums(s3[, c("coronavirus_prime_fear", "terror_prime_fear", "drunk_driving_prime_fear")], na.rm = T)
s3$priming_fear <- ifelse(s3$priming_fear==0 & s3$Group != "Control", NA, s3$priming_fear)
table(s3$priming_fear, exclude = NULL)

#### Manipulation Check ####

table(s3$mc_1, exclude = NULL)
table(s3$mc_2, exclude = NULL)

s3$Manipulation_Check <- NA
s3$Manipulation_Check[s3$mc_1==6 & s3$mc_2==8] <- 1
s3$Manipulation_Check <- ifelse(is.na(s3$Manipulation_Check), 0, s3$Manipulation_Check)
table(s3$Manipulation_Check, exclude = NULL)

s3$t1_v_control <- NA
s3$t1_v_control[s3$Group=="Treatment 1 - Coronavirus Disease"] <- "Treatment 1 - Coronavirus Disease"
s3$t1_v_control[s3$Group=="Control"] <- "Control"
table(s3$t1_v_control, exclude = NULL)

s3$t2_v_control <- NA
s3$t2_v_control[s3$Group=="Treatment 2 - Drunk Driving Crash"] <- "Treatment 2 - Drunk Driving Crash"
s3$t2_v_control[s3$Group=="Control"] <- "Control"
table(s3$t2_v_control, exclude = NULL)

s3$t3_v_control <- NA
s3$t3_v_control[s3$Group=="Treatment 3 - Terrorist Attack"] <- "Treatment 3 - Terrorist Attack"
s3$t3_v_control[s3$Group=="Control"] <- "Control"
table(s3$t3_v_control, exclude = NULL)

bartlett.test(s3$Manipulation_Check ~ s3$t1_v_control)
t.test(s3$Manipulation_Check ~ s3$t1_v_control, var.equal=T, alternative = "less")
bartlett.test(s3$Manipulation_Check ~ s3$t3_v_control)
t.test(s3$Manipulation_Check ~ s3$t3_v_control, var.equal=F, alternative = "less")
bartlett.test(s3$Manipulation_Check ~ s3$t2_v_control)
t.test(s3$Manipulation_Check ~ s3$t2_v_control, var.equal=T, alternative = "less")

table(s3$`timer_lr_placements_Page Submit`, exclude=NULL)

str(s3$`timer_treatment_1_Page Submit`)
table(s3$`timer_treatment_1_Page Submit`, exclude=NULL)
table(s3$`timer_treatment_2_Page Submit`, exclude=NULL)
table(s3$`timer_treatment_3_Page Submit`, exclude=NULL)

s3$timer_threats <- rowSums(s3[, c("timer_treatment_1_Page Submit",
                                   "timer_treatment_2_Page Submit", "timer_treatment_3_Page Submit")], na.rm = T)

table(s3$timer_threats, exclude = NULL)

table(s3$`timer_big5op_b1_Page Submit`, exclude=NULL)

s3$timer_pre_during_post_threats <- rowSums(s3[, c("timer_lr_placements_Page Submit",
                                                   "timer_threats",
                                                   "timer_big5op_b1_Page Submit")], na.rm = T)

table(s3$timer_pre_during_post_threats, exclude=NULL)

table(s3$`timer_gender_Page Submit`, exclude=NULL)
table(s3$`timer_age_Page Submit`, exclude=NULL)
table(s3$`timer_education_Page Submit`, exclude=NULL)
table(s3$`timer_income_Page Submit`, exclude=NULL)
table(s3$`timer_pid1_Page Submit`, exclude=NULL)

s3$timer_all <- rowSums(s3[, c("timer_gender_Page Submit",
                               "timer_age_Page Submit",
                               "timer_education_Page Submit",
                               "timer_income_Page Submit",
                               "timer_pid1_Page Submit",
                               "timer_pre_during_post_threats")], na.rm = T)

table(s3$timer_all, exclude = NULL)

s3$attention_bias <- (s3$timer_pre_during_post_threats/s3$timer_all)*100
table(s3$attention_bias, exclude = NULL)

s3$t1_v_control_att <- NA
s3$t1_v_control_att[s3$Group=="Treatment 1 - Coronavirus Disease"] <- "Treatment 1 - Coronavirus Disease"
s3$t1_v_control_att[s3$Group=="Control"] <- "Control"
table(s3$t1_v_control_att, exclude = NULL)

s3$t1_v_control_att <- factor(s3$t1_v_control_att, levels = c("Treatment 1 - Coronavirus Disease", "Control"), ordered = F)

s3$t2_v_control_att <- NA
s3$t2_v_control_att[s3$Group=="Treatment 2 - Drunk Driving Crash"] <- "Treatment 2 - Drunk Driving Crash"
s3$t2_v_control_att[s3$Group=="Control"] <- "Control"
table(s3$t2_v_control_att, exclude = NULL)

s3$t2_v_control_att <- factor(s3$t2_v_control_att, levels = c("Treatment 2 - Drunk Driving Crash", "Control"), ordered = F)


s3$t3_v_control_att <- NA
s3$t3_v_control_att[s3$Group=="Treatment 3 - Terrorist Attack"] <- "Treatment 3 - Terrorist Attack"
s3$t3_v_control_att[s3$Group=="Control"] <- "Control"
table(s3$t3_v_control_att, exclude = NULL)

s3$t3_v_control_att <- factor(s3$t3_v_control_att, levels = c("Treatment 3 - Terrorist Attack", "Control"), ordered = F)


bartlett.test(s3$attention_bias ~ s3$t1_v_control_att)
bartlett.test(s3$attention_bias ~ s3$t2_v_control_att)
bartlett.test(s3$attention_bias ~ s3$t3_v_control_att)

t.test(s3$attention_bias ~ s3$t1_v_control_att,  var.equal=TRUE) 
t.test(s3$attention_bias ~ s3$t3_v_control_att,  var.equal=TRUE)
t.test(s3$attention_bias ~ s3$t2_v_control_att,  var.equal=TRUE) 


#### Supplementary Appendix, Table 15. Study 3 descriptive statistics. ####

library(table1)

my.render.cont <- function(x) {
  with(stats.apply.rounding(stats.default(x), digits=2), c("",
                                                           "Mean (SD)"=sprintf("%s (&plusmn; %s)", MEAN, SD)))
}
my.render.cat <- function(x) {
  c("", sapply(stats.default(x), function(y) with(y,
                                                  sprintf("%d (%0.0f %%)", FREQ, PCT))))
}

s3$Group_Ordered <- factor(s3$Group, levels = c("Control",
                                                "Treatment 1 - Coronavirus Disease",
                                                "Treatment 3 - Terrorist Attack", 
                                                "Treatment 2 - Drunk Driving Crash"))

table1(~ Self_Conservatism + priming_fear + sex + age + education + self_left_right | Group_Ordered, data=s3, 
       render.continuous=my.render.cont, render.categorical=my.render.cat, overall=NULL)

treatment <- aov(Self_Conservatism ~ factor(Group), data = s3)
summary(treatment)
TukeyHSD(treatment)

sex <- aov(sex ~ factor(Group), data = s3)
summary(sex)

age <- aov(age ~ factor(Group), data = s3)
summary(age)

education <- aov(education ~ factor(Group), data = s3)
summary(education)

self_left_right <- aov(self_left_right ~ factor(Group), data = s3)
summary(self_left_right)
TukeyHSD(self_left_right)

#### Standardization ####

s3 <- subset(s3, subset = !is.na(Self_Conservatism))

mean(s3$Coronavirus_Prime, na.rm = T)
s3$beta_Coronavirus_Prime <- (s3$Coronavirus_Prime - mean(s3$Coronavirus_Prime, na.rm = T))
table(s3$beta_Coronavirus_Prime, exclude = NULL)

mean(s3$Terror_Prime, na.rm = T)
s3$beta_Terror_Prime <- (s3$Terror_Prime - mean(s3$Terror_Prime, na.rm = T))
table(s3$beta_Terror_Prime, exclude = NULL)

mean(s3$Drunk_Driving_Prime, na.rm = T)
s3$beta_Drunk_Driving_Prime <- (s3$Drunk_Driving_Prime - mean(s3$Drunk_Driving_Prime, na.rm = T))
table(s3$beta_Drunk_Driving_Prime, exclude = NULL)

mean(s3$sex, na.rm = T)
s3$beta_sex <- (s3$sex - mean(s3$sex, na.rm = T))
table(s3$beta_sex, exclude = NULL)

mean(s3$age, na.rm = T)
sd(s3$age, na.rm = T)
s3$beta_age <- (s3$age - mean(s3$age, na.rm = T))/(2*sd(s3$age, na.rm = T))
table(s3$beta_age, exclude = NULL)

mean(s3$education, na.rm = T)
sd(s3$education, na.rm = T)
s3$beta_education <- (s3$education - mean(s3$education, na.rm = T))/(2*sd(s3$education, na.rm = T))
table(s3$beta_education, exclude = NULL)

mean(s3$self_left_right, na.rm = T)
sd(s3$self_left_right, na.rm = T)
s3$beta_self_left_right <- (s3$self_left_right - mean(s3$self_left_right, na.rm = T))/(2*sd(s3$self_left_right, na.rm = T))
table(s3$beta_self_left_right, exclude = NULL)

mean(s3$Self_Conservatism, na.rm = T)
sd(s3$Self_Conservatism, na.rm = T)
s3$beta_Self_Conservatism <- (s3$Self_Conservatism - mean(s3$Self_Conservatism, na.rm = T))/(2*sd(s3$Self_Conservatism, na.rm = T))
table(s3$beta_Self_Conservatism, exclude = NULL)

mean(s3$priming_fear, na.rm = T)
sd(s3$priming_fear, na.rm = T)
s3$beta_priming_fear <- (s3$priming_fear - mean(s3$priming_fear, na.rm = T))/(2*sd(s3$priming_fear, na.rm = T))
table(s3$beta_priming_fear, exclude = NULL)

table(s3$State, exclude = NULL)

s3$region <- NA
s3$region[s3$State=='Istanbul'] <- "Istanbul"
s3$region[s3$State=='Tekirdag' |
            s3$State=='Edirne' |
            s3$State== 'Kirklareli' |
            s3$State== 'Balikesir' | 
            s3$State== 'Canakkale'] <- "West Marmara"
s3$region[s3$State=='Izmir' |
            s3$State== 'Aydin' | 
            s3$State==  'Denizli' |
            s3$State==  'Mugla' |
            s3$State==  'Manisa' |
            s3$State==  'Afyon' |
            s3$State== 'Kutahya' |
            s3$State== 'Usak'] <- "Aegean"
s3$region[s3$State=='Bursa' |
            s3$State== 'Eskisehir' |
            s3$State== 'Bilecik' |
            s3$State== 'Kocaeli' |
            s3$State== 'Sakarya' |
            s3$State==  'Duzce' |
            s3$State==  'Bolu' |
            s3$State== 'Yalova'] <- "East Marmara"
s3$region[s3$State=='Ankara' |
            s3$State==  'Konya' | s3$State== 'Karaman'] <- "West Anatolia"
s3$region[s3$State=='Antalya' |
            s3$State== 'Isparta' |
            s3$State==  'Adana' |
            s3$State==   'Mersin' |
            s3$State==   'Hatay' |
            s3$State==  'Kahramanmaras' |
            s3$State==   'Osmaniye' | s3$State== 'Burdur'] <- "Mediterranean"
s3$region[s3$State=='Kirikkale' |
            s3$State== 'Aksaray' |
            s3$State==  'Nigde' |
            s3$State==   'Nevsehir' |
            s3$State==   'Kirsehir' |
            s3$State==   'Kayseri' | 
            s3$State==  'Sivas' |
            s3$State==  'Yozgat'] <- "Central Anatolia"
s3$region[s3$State=='Zonguldak' |
            s3$State== 'Karabuk' |
            s3$State==  'Bartin' |
            s3$State==  'Kastamonu' |
            s3$State== 'Cankiri' |
            s3$State==  'Sinop' |
            s3$State==  'Samsun' |
            s3$State==  'Tokat' |
            s3$State==  'Corum' |
            s3$State== 'Amasya'] <- "West Black Sea"
s3$region[s3$State=='Ordu' |
            s3$State=='Trabzon' |
            s3$State=='Giresun' |
            s3$State== 'Rize' |
            s3$State== 'Artvin' |
            s3$State== 'Gumushane'] <- "East Black Sea"
s3$region[s3$State=='Erzurum' |
            s3$State=='Erzincan' |
            s3$State== 'Agri' |
            s3$State== 'Kars' | s3$State== 'Ardahan' | s3$State== 'Bayburt' | s3$State== 'Igdir'] <- "Northeast Anatolia"
s3$region[s3$State=='Malatya' |
            s3$State== 'Elazig' |
            s3$State== 'Bingol' |
            s3$State== 'Tunceli' |
            s3$State== 'Van' |
            s3$State== 'Mus' |
            s3$State== 'Bitlis' |
            s3$State== 'Hakkari'] <- "Central East Anatolia"
s3$region[s3$State=='Gaziantep' |
            s3$State== 'Adiyaman' |
            s3$State== 'Kilis' |
            s3$State== 'Sanliurfa' |
            s3$State== 'Diyarbakir' |
            s3$State== 'Mardin' |
            s3$State== 'Batman' |
            s3$State== 'Sirnak' |
            s3$State== 'Siirt'] <- "Southeast Anatolia"
table(s3$region, exclude = NULL)
table(s3$State, s3$region, exclude = NULL)
s3$region <- ifelse(is.na(s3$region), "Outside Turkey", s3$region)

#### Supplementary Appendix, Table 18. Study 3 balance testing. ####

library(nnet)

balance_test <-  multinom(factor(Group) ~ beta_sex + beta_age + beta_education + beta_self_left_right + beta_akp + beta_chp + beta_other_right + beta_other_left, data=s3)
summary(balance_test)
lmtest::lrtest(balance_test)
tab_model(balance_test, show.ci = 0.95)

stargazer(balance_test, type="html", out="s3_balance_test.htm",
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          covariate.labels =c("Sex: Male", "Age", "Education", "Left-Right Ideological Placement", "Partisanship: AKP", "Partisanship: CHP", "Partisanship: Other Parties on the Right", "Partisanship: Other Parties on the Left", "Intercept"),
          dep.var.labels   = c("Coronavirus MS", "Terrorism MS", "Drunk Driving Crash MS"))

#### Supplementary Appendix, Table 22. Study 3 causal mediation analysis models. ####

#### Model 1 ####

out.fit.m01 <- lm(beta_Self_Conservatism ~ beta_Coronavirus_Prime + beta_Terror_Prime + beta_Drunk_Driving_Prime, data=s3)
summary(out.fit.m01)
tab_model(out.fit.m01, show.ci = 0.95)

confint(out.fit.m01, 'beta_Coronavirus_Prime', level=0.95)
confint(out.fit.m01, 'beta_Terror_Prime', level=0.95)
confint(out.fit.m01, 'beta_Drunk_Driving_Prime', level=0.95)

C <- diag(4)
model.mc <- glht(out.fit.m01, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

med.fit.m1 <- lm(beta_priming_fear ~ beta_Coronavirus_Prime + beta_Terror_Prime + beta_Drunk_Driving_Prime, data=s3)
summary(med.fit.m1)
tab_model(med.fit.m1, show.ci = 0.95)

C <- diag(4)
model.mc <- glht(med.fit.m1, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

out.fit.m1 <- lm(beta_Self_Conservatism ~ beta_Coronavirus_Prime + beta_Terror_Prime + beta_Drunk_Driving_Prime + beta_priming_fear, data=s3)
summary(out.fit.m1)
tab_model(out.fit.m1, show.ci = 0.95)

confint(out.fit.m1, 'beta_Coronavirus_Prime', level=0.95)
confint(out.fit.m1, 'beta_Terror_Prime', level=0.95)
confint(out.fit.m1, 'beta_Drunk_Driving_Prime', level=0.95)

C <- diag(5)
model.mc <- glht(out.fit.m1, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

library(mediation)

set.seed(123)
med.out.m1.coronavirus <- mediate(med.fit.m1, out.fit.m1, treat = "beta_Coronavirus_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.out.m1.coronavirus)
med.out.m1.coronavirus$d.avg.ci

psych::describe(s3$beta_Coronavirus_Prime)
psych::describe(s3$beta_priming_fear)
psych::describe(s3$beta_Self_Conservatism)

set.seed(123)
med.out.m1.terror <- mediate(med.fit.m1, out.fit.m1, treat = "beta_Terror_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.out.m1.terror)
med.out.m1.terror$d.avg.ci

psych::describe(s3$beta_Terror_Prime)
psych::describe(s3$beta_priming_fear)
psych::describe(s3$beta_Self_Conservatism)

set.seed(123)
med.out.m1.general <- mediate(med.fit.m1, out.fit.m1, treat = "beta_Drunk_Driving_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.out.m1.general)
med.out.m1.general$d.avg.ci

psych::describe(s3$beta_Drunk_Driving_Prime)
psych::describe(s3$beta_priming_fear)
psych::describe(s3$beta_Self_Conservatism)

#### Model 2 ####

out.fit.m02 <- lm(beta_Self_Conservatism ~ beta_Coronavirus_Prime + beta_Terror_Prime + beta_Drunk_Driving_Prime + beta_sex + beta_age + beta_education + beta_self_left_right, data=s3)
summary(out.fit.m02)
tab_model(out.fit.m02, show.ci = 0.95)

confint(out.fit.m02, 'beta_Coronavirus_Prime', level=0.95)
confint(out.fit.m02, 'beta_Terror_Prime', level=0.95)
confint(out.fit.m02, 'beta_Drunk_Driving_Prime', level=0.95)

C <- diag(8)
model.mc <- glht(out.fit.m02, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

med.fit.m2 <- lm(beta_priming_fear ~ beta_Coronavirus_Prime + beta_Terror_Prime + beta_Drunk_Driving_Prime + beta_sex + beta_age + beta_education + beta_self_left_right, data=s3)
summary(med.fit.m2)
tab_model(med.fit.m2, show.ci = 0.95)

C <- diag(8)
model.mc <- glht(med.fit.m2, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

out.fit.m2 <- lm(beta_Self_Conservatism ~ beta_Coronavirus_Prime + beta_Terror_Prime + beta_Drunk_Driving_Prime + beta_priming_fear + beta_sex + beta_age + beta_education + beta_self_left_right, data=s3)
summary(out.fit.m2)
tab_model(out.fit.m2, show.ci = 0.95)

confint(out.fit.m2, 'beta_Coronavirus_Prime', level=0.95)
confint(out.fit.m2, 'beta_Terror_Prime', level=0.95)
confint(out.fit.m2, 'beta_Drunk_Driving_Prime', level=0.95)

C <- diag(9)
model.mc <- glht(out.fit.m2, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

library(mediation)

set.seed(123)
med.out.m2.coronavirus <- mediate(med.fit.m2, out.fit.m2, treat = "beta_Coronavirus_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.out.m2.coronavirus)

set.seed(123)
med.out.m2.terror <- mediate(med.fit.m2, out.fit.m2, treat = "beta_Terror_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.out.m2.terror)

set.seed(123)
med.out.m2.general <- mediate(med.fit.m2, out.fit.m2, treat = "beta_Drunk_Driving_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.out.m2.general)

library(stargazer)

stargazer(out.fit.m01, med.fit.m1, out.fit.m1, out.fit.m02, med.fit.m2, out.fit.m2,
          type = "html", title=" ", digits=2, out="s3_model_outputs.htm",
          star.cutoffs = c(0.05, 0.01, 0.001),
          model.numbers = F,
          column.labels = c("Model 0.1.1", "Model 1.1", "Model 1.2", "Model 0.2.1", "Model 2.1", "Model 2.2"),
          covariate.labels = c("Coronavirus Mortality Salience", 
                               "Terrorism Mortality Salience", 
                               "Drunk Driving Crash Mortality Salience", 
                               "Mortality Fear",
                               "Sex: Male",
                               "Age", 
                               "Education", 
                               "Left-Right Ideological Placement"),
          add.lines = list(c("Estimation", "OLS", "OLS", "OLS", "OLS", "OLS", "OLS")))

#### Unreported Moderated Mediation by Sex ####
med.fit.m1.sex <- lm(beta_priming_fear ~ beta_sex*beta_Coronavirus_Prime + beta_sex*beta_Terror_Prime + beta_sex*beta_Drunk_Driving_Prime, data=s3)
summary(med.fit.m1.sex)
tab_model(med.fit.m1.sex, show.ci = 0.95)

C <- diag(8)
model.mc <- glht(med.fit.m1.sex, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

out.fit.m1.sex <- lm(beta_Self_Conservatism ~ beta_sex*beta_Coronavirus_Prime + beta_sex*beta_Terror_Prime + beta_sex*beta_Drunk_Driving_Prime + beta_sex*beta_priming_fear, data=s3)
summary(out.fit.m1.sex)
tab_model(out.fit.m1.sex, show.ci = 0.95)

C <- diag(10)
model.mc <- glht(out.fit.m1.sex, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
round(BH.OUT$test$pvalues, 3)

library(mediation)

table(s3$beta_sex)

set.seed(123)
med.fit.m1.coronavirus.sex <- mediate(med.fit.m1.sex, out.fit.m1.sex, treat = "beta_Coronavirus_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.fit.m1.coronavirus.sex)
test.modmed(med.fit.m1.coronavirus.sex, covariates.1 = list(beta_sex = -0.558504221954162), covariates.2 = list(beta_sex = 0.441495778045838), sims = 10000)

set.seed(123)
med.fit.m1.terror.sex <- mediate(med.fit.m1.sex, out.fit.m1.sex, treat = "beta_Terror_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.fit.m1.terror.sex)
test.modmed(med.fit.m1.terror.sex, covariates.1 = list(beta_sex = -0.558504221954162), covariates.2 = list(beta_sex = 0.441495778045838), sims = 10000)

set.seed(123)
med.fit.m1.general.sex <- mediate(med.fit.m1.sex, out.fit.m1.sex, treat = "beta_Drunk_Driving_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.fit.m1.general.sex)
test.modmed(med.fit.m1.general.sex, covariates.1 = list(beta_sex = -0.558504221954162), covariates.2 = list(beta_sex = 0.441495778045838), sims = 10000)

#### Unreported Moderated Mediation by Age ####
med.fit.m1.age <- lm(beta_priming_fear ~ beta_age*beta_Coronavirus_Prime + beta_age*beta_Terror_Prime + beta_age*beta_Drunk_Driving_Prime, data=s3)
summary(med.fit.m1.age)
tab_model(med.fit.m1.age, show.ci = 0.95)

C <- diag(8)
model.mc <- glht(med.fit.m1.age, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

out.fit.m1.age <- lm(beta_Self_Conservatism ~ beta_age*beta_Coronavirus_Prime + beta_age*beta_Terror_Prime + beta_age*beta_Drunk_Driving_Prime + beta_age*beta_priming_fear, data=s3)
summary(out.fit.m1.age)
tab_model(out.fit.m1.age, show.ci = 0.95)

C <- diag(10)
model.mc <- glht(out.fit.m1.age, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

table(s3$beta_age, exclude = NULL)

library(mediation)

set.seed(123)
med.fit.m1.coronavirus.age <- mediate(med.fit.m1.age, out.fit.m1.age, treat = "beta_Coronavirus_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.fit.m1.coronavirus.age)
test.modmed(med.fit.m1.coronavirus.age, covariates.1 = list(beta_age = -1.32755860717187), covariates.2 = list(beta_age = 1.81179161707738), sims = 10000)

set.seed(123)
med.fit.m1.terror.age <- mediate(med.fit.m1.age, out.fit.m1.age, treat = "beta_Terror_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.fit.m1.terror.age)
test.modmed(med.fit.m1.terror.age, covariates.1 = list(beta_age = -1.32755860717187), covariates.2 = list(beta_age = 1.81179161707738), sims = 10000)

set.seed(123)
med.fit.m1.general.age <- mediate(med.fit.m1.age, out.fit.m1.age, treat = "beta_Drunk_Driving_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.fit.m1.general.age)
test.modmed(med.fit.m1.general.age, covariates.1 = list(beta_age = -1.32755860717187), covariates.2 = list(beta_age = 1.81179161707738), sims = 10000)

#### Unreported Moderated Mediation by Left-Right Ideological Placement ####
med.fit.m1.ideology <- lm(beta_priming_fear ~ beta_self_left_right*beta_Coronavirus_Prime + beta_self_left_right*beta_Terror_Prime + beta_self_left_right*beta_Drunk_Driving_Prime, data=s3)
summary(med.fit.m1.ideology)
tab_model(med.fit.m1.ideology, show.ci = 0.95)

C <- diag(8)
model.mc <- glht(med.fit.m1.ideology, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

out.fit.m1.ideology <- lm(beta_Self_Conservatism ~ beta_self_left_right*beta_Coronavirus_Prime + beta_self_left_right*beta_Terror_Prime + beta_self_left_right*beta_Drunk_Driving_Prime + beta_self_left_right*beta_priming_fear, data=s3)
summary(out.fit.m1.ideology)
tab_model(out.fit.m1.ideology, show.ci = 0.95)

C <- diag(10)
model.mc <- glht(out.fit.m1.ideology, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

table(s3$beta_self_left_right, exclude = NULL)

library(mediation)

set.seed(123)
med.fit.m1.coronavirus.ideology <- mediate(med.fit.m1.ideology, out.fit.m1.ideology, treat = "beta_Coronavirus_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.fit.m1.coronavirus.ideology)
test.modmed(med.fit.m1.coronavirus.ideology, covariates.1 = list(beta_self_left_right = -0.64007101121113), covariates.2 = list(beta_self_left_right = 1.04070225539218), sims = 10000)

set.seed(123)
med.fit.m1.terror.ideology <- mediate(med.fit.m1.ideology, out.fit.m1.ideology, treat = "beta_Terror_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.fit.m1.terror.ideology)
test.modmed(med.fit.m1.terror.ideology, covariates.1 = list(beta_self_left_right = -0.64007101121113), covariates.2 = list(beta_self_left_right = 1.04070225539218), sims = 10000)

set.seed(123)
med.fit.m1.general.ideology <- mediate(med.fit.m1.ideology, out.fit.m1.ideology, treat = "beta_Drunk_Driving_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.fit.m1.general.ideology)
test.modmed(med.fit.m1.general.ideology, covariates.1 = list(beta_self_left_right = -0.64007101121113), covariates.2 = list(beta_self_left_right = 1.04070225539218), sims = 10000)

#### Unreported Serial Mediation Analyses (provided for the interested readers) ####
##### Note to the Readers #####

# Serial mediation analyses, testing the alternative causal chain running from treatments to mortality fear to personality traits to conservatism, 
# revealed that the theorized mechanisms were confounded if I account for the alternative paths running through the personality trait Openness 
# However, there was also no significant and replicable mediating effect operating through Openness. 
# I did not consider the experimentally induced mortality fear’s confounded mediating effect in the serial mediation analysis to be 
# of any theoretical significance to discuss further. This is because the personality trait Openness remained stable in the design, and 
# I did not theorize any implications of personality traits for the causal mechanism as mediators.

names(s3)[names(s3)=="big5op_1"] <- "openness_1"
names(s3)[names(s3)=="big5op_2"] <- "openness_2"
names(s3)[names(s3)=="big5op_3"] <- "openness_3"
names(s3)[names(s3)=="big5op_4"] <- "openness_4"
names(s3)[names(s3)=="big5op_5"] <- "openness_5"
names(s3)[names(s3)=="big5op_6"] <- "openness_6"
s3$big5op_7 <- 8 - s3$big5op_7
names(s3)[names(s3)=="big5op_7"] <- "openness_7"
names(s3)[names(s3)=="big5op_8"] <- "openness_8"
s3$big5op_9 <- 8 - s3$big5op_9
names(s3)[names(s3)=="big5op_9"] <- "openness_9"
names(s3)[names(s3)=="big5op_10"] <- "openness_10"

psych::alpha(s3[,c("openness_1", "openness_2", "openness_3", "openness_4",
                   "openness_5", "openness_6", "openness_7", "openness_8", "openness_9", "openness_10")])

s3$openness <- s3$openness_1 + 
  s3$openness_2 + 
  s3$openness_3 + 
  s3$openness_4 + 
  s3$openness_5 + 
  s3$openness_6 + 
  s3$openness_7 + 
  s3$openness_8 +
  s3$openness_9 +
  s3$openness_10

psych::describeBy(s3$openness, s3$Group)

mean(s3$openness, na.rm = T)
sd(s3$openness, na.rm = T)
s3$beta_openness <- (s3$openness - mean(s3$openness, na.rm = T))/(2*sd(s3$openness, na.rm = T))
table(s3$beta_openness, exclude = NULL)

library(lavaan)

model=
  "
  #Regressions
  beta_priming_fear ~ tc_f*beta_Coronavirus_Prime + tt_f*beta_Terror_Prime + td_f*beta_Drunk_Driving_Prime + beta_sex + beta_age + beta_education + beta_self_left_right
  beta_openness ~ f_o*beta_priming_fear + tc_o*beta_Coronavirus_Prime + tt_o*beta_Terror_Prime + td_o*beta_Drunk_Driving_Prime + beta_sex + beta_age + beta_education + beta_self_left_right
  beta_Self_Conservatism ~ op*beta_openness + f*beta_priming_fear + c*beta_Coronavirus_Prime + t*beta_Terror_Prime + d*beta_Drunk_Driving_Prime + beta_sex + beta_age + beta_education + beta_self_left_right
  
  #Defined Parameters:
  
  ie_all_c_o := tc_f*f_o*op
  ie_all_t_o := tt_f*f_o*op
  ie_all_d_o := td_f*f_o*op
  
  
  ie_short_c_o := tc_o*op
  ie_short_t_o := tt_o*op
  ie_short_d_o := td_o*op
  
  ie_c_f := tc_f*f
  ie_t_f := tt_f*f
  ie_d_f := td_f*f
  
  de_c := c
  de_t := t
  de_d := d

"

fit <-  sem(model, data=s3, se = "bootstrap", bootstrap = 10000)
summary(fit, fit.measures=TRUE, standardized=TRUE, rsquare=TRUE, nd = 2)
parameterEstimates(fit, rsquare=T, output = "text", level = 0.90)

psych::describe(s3$beta_Terror_Prime)
psych::describe(s3$beta_priming_fear)
psych::describe(s3$beta_openness)
psych::describe(s3$beta_Self_Conservatism)

#### Unreported Study 3 descriptive statistics of the personality trait Openness. ####

table1(~ openness | Group_Ordered, data=s3, 
       render.continuous=my.render.cont, render.categorical=my.render.cat, overall=NULL)

#### Supplementary Appendix, Study 3 Sample Representativeness Tests, Tables 9 to 12. ####

table(s3$age, exclude = NULL)
s3$age_cat <- NA
s3$age_cat[s3$age>=18 & s3$age<=24] <- "18-24"
s3$age_cat[s3$age>=25 & s3$age<=34] <- "25-34"
s3$age_cat[s3$age>=35 & s3$age<=49] <- "35-49"
s3$age_cat[s3$age>=50 & s3$age<=64] <- "50-64"
s3$age_cat[s3$age>=65] <- "65 or Over"
table(s3$age_cat, exclude = NULL)
round(prop.table(table(s3$age_cat, exclude = NULL))*100, 1)

s3_age <- c(4.8, 6.1, 20.5, 53.1, 15.4)
tr_2019_age <- c(15.2, 21, 30.1, 21.1, 12.6)

age_s3_tr_2018 <- cbind(s3_age, tr_2019_age)
row.names(age_s3_tr_2018) <- c("18-24","25-34", "35-49", "50-64","65+")

chisq.test(age_s3_tr_2018)
round(chisq.test(age_s3_tr_2018 )$stdres, 2)
round(rcompanion::cramerV(age_s3_tr_2018), 2)

table(s3$sex, exclude = NULL)
s3$sex_cat <- NA
s3$sex_cat[s3$sex==1] <- "Male"
s3$sex_cat[s3$sex==0] <- "Female"
table(s3$sex_cat, exclude = NULL)
round(prop.table(table(s3$sex_cat, exclude = NULL))*100, 1)

s3_gender <- c(55.8, #male
               44.1) #female

tr_2019_gender <- c(50.7,
                    49.3)

gender_s3_tr_2018 <- cbind(s3_gender, tr_2019_gender)
row.names(gender_s3_tr_2018) <- c("Male","Female")

chisq.test(gender_s3_tr_2018)
round(rcompanion::cramerV(gender_s3_tr_2018), 2)

table(s3$education, exclude = NULL)
s3$education_cat <- NA
s3$education_cat[s3$education==1 | s3$education==2 | s3$education==3] <- "Below High School Graduate"
s3$education_cat[s3$education==4 | s3$education==5] <- "High School Graduate"
s3$education_cat[s3$education==6 | s3$education==7 | s3$education==8] <- "College Graduate or Above"
s3$education_cat[is.na(s3$education)] <- "Unknown"
table(s3$education_cat, exclude = NULL)
s3$education_cat <- factor(s3$education_cat, levels = c("Below High School Graduate", "High School Graduate",
                                                        "College Graduate or Above", "Unknown"))
round(prop.table(table(s3$education_cat, exclude = NULL))*100, 1)

s3_education <- c(11.6, 31.1, 57.3)
tr_2019_education <- c(46.8, 25.1, 18.9)

education_s3_tr_2018 <- cbind(s3_education, tr_2019_education)
row.names(education_s3_tr_2018) <- c("Below High School Graduate",
                                     "High School Graduate",
                                     "College Graduate or Above")

chisq.test(education_s3_tr_2018)
round(chisq.test(education_s3_tr_2018 )$stdres, 2)
round(rcompanion::cramerV(education_s3_tr_2018), 2)

round(prop.table(table(s3$State, exclude = NULL))*100, 1)

s3_geo <- c(0.7,
            0.8,
            0.1,
            0.4,
            1,
            1.9,
            1.7,
            12,
            0.6,
            0.6,
            0,
            0.1,
            0.1,
            0.2,
            0.1,
            0,
            0.2,
            0,
            0,
            0.1,
            0,
            0,
            0,
            0,
            0.4,
            0.1,
            0.1,
            0.4,
            0,
            0,
            0,
            1.6,
            0.4,
            1,
            0.1,
            1.6,
            5.9,
            0.1,
            0.2,
            30,
            0.2,
            1.8,
            1.3,
            2.4,
            3.7,
            0,
            0,
            0.4,
            0.4,
            0.1,
            0.1,
            0.1,
            0,
            0,
            0,
            0.6,
            0.5,
            0.6,
            0.1,
            0,
            0,
            0.4,
            0.1,
            0,
            2,
            0.1,
            15.8,
            0.7,
            0.2,
            0.1,
            2,
            0.2,
            0.1,
            0.1,
            0,
            0,
            0.6,
            1.1,
            0.5,
            0.1,
            0)

tr_2019_geo <- c(1.4,
                 1.3,
                 0.5,
                 0.7,
                 0.9,
                 1.8,
                 1.3,
                 5.7,
                 1.7,
                 0.5,
                 0.4,
                 0.3,
                 0.4,
                 0.8,
                 0.5,
                 0.4,
                 1,
                 0.3,
                 0.4,
                 0.7,
                 0.3,
                 0.4,
                 0.1,
                 1.1,
                 1,
                 1,
                 0.6,
                 0.5,
                 0.2,
                 0.3,
                 0.2,
                 2.4,
                 1.3,
                 1.1,
                 0.3,
                 0.5,
                 3.8,
                 0.4,
                 0.3,
                 18.5,
                 1.3,
                 2.2,
                 1.9,
                 2.7,
                 3,
                 0.3,
                 0.6,
                 0.6,
                 0.9,
                 0.3,
                 0.3,
                 0.5,
                 0.1,
                 0.1,
                 0.2,
                 2.2,
                 1.9,
                 1.8,
                 0.7,
                 0.6,
                 0.2,
                 0.8,
                 0.3,
                 0.5,
                 2.7,
                 0.3,
                 6.9,
                 0.8,
                 0.3,
                 0.3,
                 1.7,
                 0.4,
                 0.7,
                 0.5,
                 0.3,
                 0.8,
                 1.7,
                 1.3,
                 0.7,
                 0.5,
                 0.5)

geo_s3_tr_2018 <- cbind(s3_geo, tr_2019_geo)
row.names(geo_s3_tr_2018) <- c(
  "Aydin",
  "Denizli",
  "Usak",
  "Kutahya",
  "Afyon",
  "Manisa",
  "Mugla",
  "Izmir",
  "Kayseri",
  "Aksaray",
  "Kirikkale",
  "Kirsehir",
  "Nevsehir",
  "Sivas",
  "Yozgat",
  "Nigde",
  "Malatya",
  "Bingol",
  "Bitlis",
  "Elazig",
  "Hakkari",
  "Mus",
  "Tunceli",
  "Van",
  "Ordu",
  "Trabzon",
  "Giresun",
  "Rize",
  "Artvin",
  "Bartin",
  "Gumushane",
  "Kocaeli",
  "Sakarya",
  "Eskisehir",
  "Bilecik",
  "Duzce",
  "Bursa",
  "Bolu",
  "Yalova",
  "Istanbul",
  "Kahramanmaras",
  "Mersin",
  "Hatay",
  "Adana",
  "Antalya",
  "Burdur",
  "Isparta",
  "Osmaniye",
  "Erzurum",
  "Erzincan",
  "Kars",
  "Agri",
  "Ardahan",
  "Bayburt",
  "Igdir",
  "Gaziantep",
  "Sanliurfa",
  "Diyarbakir",
  "Adiyaman",
  "Batman",
  "Kilis",
  "Mardin",
  "Siirt",
  "Sirnak",
  "Konya",
  "Karaman",
  "Ankara",
  "Zonguldak",
  "Cankiri",
  "Karabuk",
  "Samsun",
  "Amasya",
  "Corum",
  "Kastamonu",
  "Sinop",
  "Tokat",
  "Balikesir",
  "Tekirdag",
  "Canakkale",
  "Edirne",
  "Kirklareli")

chisq.test(geo_s3_tr_2018)
round(chisq.test(geo_s3_tr_2018)$stdres, 2)
round(rcompanion::cramerV(geo_s3_tr_2018), 2)

#### Sequential G-Estimation ####
##### Note to the Readers #####
# If you have already run the codes for the serial mediation analyses above, you do not need to recreate the personality trait variable Openness (lines between 1118 and 1143) below.
# Instead, please start running the codes starting from the line 1147.

names(s3)[names(s3)=="big5op_1"] <- "openness_1"
names(s3)[names(s3)=="big5op_2"] <- "openness_2"
names(s3)[names(s3)=="big5op_3"] <- "openness_3"
names(s3)[names(s3)=="big5op_4"] <- "openness_4"
names(s3)[names(s3)=="big5op_5"] <- "openness_5"
names(s3)[names(s3)=="big5op_6"] <- "openness_6"
s3$big5op_7 <- 8 - s3$big5op_7
names(s3)[names(s3)=="big5op_7"] <- "openness_7"
names(s3)[names(s3)=="big5op_8"] <- "openness_8"
s3$big5op_9 <- 8 - s3$big5op_9
names(s3)[names(s3)=="big5op_9"] <- "openness_9"
names(s3)[names(s3)=="big5op_10"] <- "openness_10"

psych::alpha(s3[,c("openness_1", "openness_2", "openness_3", "openness_4",
                   "openness_5", "openness_6", "openness_7", "openness_8", "openness_9", "openness_10")])

s3$openness <- s3$openness_1 + 
  s3$openness_2 + 
  s3$openness_3 + 
  s3$openness_4 + 
  s3$openness_5 + 
  s3$openness_6 + 
  s3$openness_7 + 
  s3$openness_8 +
  s3$openness_9 +
  s3$openness_10

psych::describeBy(s3$openness, s3$Group)

mean(s3$openness, na.rm = T)
sd(s3$openness, na.rm = T)
s3$beta_openness <- (s3$openness - mean(s3$openness, na.rm = T))/(2*sd(s3$openness, na.rm = T))

set.seed(123)

sqm1 <- lm(beta_Self_Conservatism ~ beta_Coronavirus_Prime + beta_Terror_Prime + beta_Drunk_Driving_Prime + beta_priming_fear + beta_openness, data=s3)
summary(sqm1)

s3$beta_Self_Conservatism_Adjusted <- (s3$beta_Self_Conservatism - s3$beta_priming_fear*sqm1$coefficients["beta_priming_fear"])

sge_function <- function(data, sample) {
  
  d <- data[sample,]
  
  sgm2 <- lm(beta_Self_Conservatism_Adjusted ~ beta_Coronavirus_Prime + beta_Terror_Prime + beta_Drunk_Driving_Prime + beta_sex + beta_age + beta_education + beta_self_left_right, data=d)
  summary(sgm2)
  
  sgm3 <- lm(beta_Self_Conservatism ~ beta_Coronavirus_Prime + beta_Terror_Prime + beta_Drunk_Driving_Prime + beta_sex + beta_age + beta_education + beta_self_left_right, data=d)
  summary(sgm3)
  
  return(c((sgm3$coefficients["beta_Coronavirus_Prime"] - sgm2$coefficients["beta_Coronavirus_Prime"]),
         (sgm3$coefficients["beta_Terror_Prime"] - sgm2$coefficients["beta_Terror_Prime"]),
  (sgm3$coefficients["beta_Drunk_Driving_Prime"] - sgm2$coefficients["beta_Drunk_Driving_Prime"]))) 
  
}

sgm_out <- boot::boot(data=s3, statistic=sge_function, R=10000)
sgm_out
sgm_out_coronavirus_ci <- boot::boot.ci(sgm_out, type = "bca", index = 1) # 95%   (-0.0781, -0.0685 )
sgm_out_coronavirus_ci
sgm_out_coronavirus_ci$t0 # -0.073 [-0.078, -0.068]
sgm_out_coronavirus_ci$bca
sgm_out_terror_ci <- boot::boot.ci(sgm_out, type = "bca", index = 2) # 95%   (-0.0740, -0.0634 )
sgm_out_terror_ci$t0 # -0.069 [-0.074, -0.063]
sgm_out_terror_ci$bca
sgm_out_drunk_driving_ci <- boot::boot.ci(sgm_out, type = "bca", index = 3) # 95%   (-0.0973, -0.0864 )
sgm_out_drunk_driving_ci$t0 # -0.092 [-0.097, -0.086]
sgm_out_drunk_driving_ci$bca

